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Creative telescoping algorithms compute linear differential 
equations satisfied by multiple integrals with parameters. 
We describe a precise and elementary algorithmic version 
of the GrifEths-Dwork method for the creative telescoping 
of rational functions. This leads to bounds on the order and 
degree of the coefficients of the differential equation, and to 
the first complexity result which is simply exponential in the 
number of variables. One of the important features of the al- 
gorithm is that it does not need to compute certificates. The 
approach is vindicated by a prototype implementation. 

Categories and Subject Descriptors: 
1.1.2 [Computing Methodologies]: Symbolic and Alge- 
braic Manipulations — Algebraic Algorithms 

General Terms: Algorithms, Theory. 

Keywords: Integration, creative telescoping, algorithms, 
complexity, Picard-Fuchs equation. 

1. INTRODUCTION 

In computer algebra, creative telescoping is an approach in- 
troduced by Zeilberger to address definite summation and 
integration of a large class of functions and sequences [28,29, 
27]. Its vast scope includes the computation of differential 
equations for multiple integrals with parameters of rational 
or algebraic functions. Within this class, creative telescop- 
ing is similar to well-studied older approaches whose key no- 
tion is the Picard-Fuchs differential equation, see e.g. [23]. 

We study the multivariate rational case: Given a rational 
function F(t, xi, . . . , x n ), we aim at finding n other rational 
functions Ai(t, xi, . . . , x n ) and a differential operator with 
polynomial coefficients ^ r =o c i(t)d 3 t , denoted T, such that 



T(F)^5>(t)#F = 



d t A t 
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where d{ denotes ~ and di denotes The integer r is 
the order of T and maxj deg Cj is its degree. 

Throughout the article, the constant field k of F is as- 
sumed to be of characteristic zero. Under suitable hypothe- 
ses, the operator T is a differential equation satisfied by inte- 
grals J Fdx over a domain 7, without boundaries, where F 
has no pole. A misbehavior can happen when the certifi- 
cate has poles outside those of F: it may not be possible 
to integrate term by term the right-hand side of Eq. (1), 
see §4.1. The certificate is called regular when it does not 
contain poles other than those of F. For integration, there 
is no need to compute the actual certificate provided that it 
is regular. 

Several methods are known that can find a telescoper and 
the corresponding certificate [17,26,7,15]. However, the 
practical cost of using these methods in multivariate prob- 
lems remains high and a better understanding of the size or 
complexity of the objects of creative telescoping is clearly 
needed. The present work is part of the on-going effort in 
this direction [2,3,4]. The study of the rational case is moti- 
vated both by its fundamental nature and by its applications 
to the computation of diagonals in combinatorics, number 
theory and physics [17,6,20]. The rational case with n vari- 
ables also includes the algebraic case with n — 1 variables [4]. 

Previous works. An obviously related problem is, given a ra- 
tional function F(x\, . . . ,x n ), to decide whether there exist 
rational functions A\ , . . . , A n such that F equals ^"-1 Tla^- 
When n — 1, the latter question is easily solved by Hermite 
reduction. This is the basis of an algorithm for creative 
telescoping [3] that we outline in §2.1. 

Picard [25, chap. 7] gave methods when n — 2 from which 
he deduced that a telescoping equation can always be found 
in that case [24]. This too, has led to an algorithm [4]. The 
Griffiths-Dwork method [8, §3; 9, §8; 12] solves the prob- 
lem for a general n, in the setting of de Rham cohomology 
and under a regularity assumption. The method can be 
viewed as a generalization of Hermite reduction. Indepen- 
dently, Christol used a similar method to prove that diago- 
nals of rational functions, under a regularity hypothesis, are 
differentially finite [5]; then he noticed that a deformation 
technique, for which he credits Dwork, can handle singu- 
lar cases [6]. The Griffiths-Dwork method is also used in 
point counting applications [1, 11] and the study of mirror 
maps [20]. 

In terms of complexity, not much is known. If a rational 
function F(t, xi, . . . , x„) has degree d, a study of Lipshitz's 
argument [17] shows that there exists a telescoper of order 
and degree with a regular certificate of size d 0( - n \ 



Most algorithms [17, 28, 26, 7, 2, 15] cannot avoid the com- 
putation of the certificate, which impacts their complexity. 
The complexity of Lipshitz's algorithm is d°^ n ' operations 
in k; the complexity of no other algorithm is known. Pan- 
cratz [22] developed an approach similar to ours, under a re- 
strictive hypothesis, much stronger than Griffiths' regularity 
assumption. He proceeds to a complexity analysis of his al- 
gorithm but in terms of operations in k(t) rather than in the 
base field k. Algorithms based on non-commutative Grob- 
ner bases and elimination [28, 26] or based on the search of 
rational solutions to differential equations [7] resist to com- 
plexity analysis. The method of Apagodu and Zeilberger [2] 
requires a generic exponent and specialization seems prob- 
lematic. 

For the restricted class of diagonals of rational functions, 
there is a heuristic based on series expansion and differential 
approximation [14] ; it does not need to compute a certificate. 
However, even using the bounds m^d°^- n \ its direct imple- 
mentation has a complexity of d° (n - 1 operations in k. 

Contributions. Our main result, obtained with the Grifliths- 
Dwork method and a deformation technique, is the existence 
of a telescoper with regular certificate of order at most d n 
and degree d°^ that can be computed in d°^ arithmetic 
operations in k. For generic rational functions, the tele- 
scoper computed is the minimal order telescoper with regu- 
lar certificate. Theorems 6, 10 and 11 state precise estimates. 
To the best of our knowledge, the bounds on the order and 
degree are better than what was known and it is the first 
time that a complexity simply exponential in n is reached. 
Note that we do not know whether a rational function ad- 
mits a telescoper whose certificate has size polynomial in d n , 
but our algorithm does not need to compute the certificate. 
A prototype implementation shows that this algorithm can 
lead to a spectacular improvement over previous methods, 
though the domain of improvement is not satisfactory yet. 

Acknowledgement. We are deeply grateful to Gilles Christol 
for many rewarding discussions. 

2. OVERVIEW OF THE METHOD 

In this section we introduce the basics of the Griffiths-Dwork 
method. In dimension 1, this method coincides with classical 
Hermite reduction which we first recall. 

2.1 Dimension one: Hermite reduction 

Let F be a rational function in x, over a field L, written 
as a/f e , with a and / two polynomials, the latter being 
square- free, i.e. the polynomials d x f and / are coprime. In 
particular a equals uf + vd x f for some polynomials u and v. 
Then, if i > 1, the function F rewrites 

U + j^-rdxV 
F = f \ +3. 



(t-l)f* 



Iterating this reduction step £ times gives F as y + d x jrrr 
for some polynomials U and V. Next, Euclidean division 
allows to write U as r + sf, with r of degree less than the 
degree of /, yielding the additive decomposition 



F 



f 



V 



f 



denoted [F]. This reduced form features important proper- 
ties: 

(Linearity) f being fixed, [F] depends linearly on F; 

(Soundness) if [F] is zero, then F is a derivative w.r.t. x; 

(Confinement) [F] lies in a finite-dimensional vector space 
over L depending only on / (with dimension deg^, /); 

(Normalization) if F is a derivative w.r.t. x, then [F] is zero. 

These properties are enough to compute a telescoper: As- 
sume now that L is k(t) for a field k. If for some elements 
of L, say ao, . . . , a p , the reduced form [y^ aidlF^ vanishes 
then the operator y\ aid\ is a telescoper, thanks to the 
soundness property. Thanks to the linearity property, this 
is equivalent to the vanishing of y\ a; \d\F~\ . Thanks to 
the confinement property, it is always possible to find such 
a relation. Thanks to the normalization property, every tele- 
scoper arises this way. In particular, so does the telescoper 
of minimal order. 

2.2 Higher dimension: Griffiths-Dwork 
reduction 

Let F be a rational function in n variables Xx, . . . , x n , writ- 
ten as a/f , with / a square- free polynomial. If I > 1 and 
if a lies in the ideal generated by / and its derivatives dif , 
then we can write a as wf + yV Vidif, for some polynomials 
it, vi, . . . , v n , and F rewrites 



f 



d, 



Provided that this ideal contains all polynomials, any F can 
be reduced to a function with simple poles by iteration of 
this identity. The soundness and linearity properties are 
naturally satisfied, but extending further the reduction to 
obtain at least the confinement property is not straightfor- 
ward and requires stronger assumptions [21]. A difficulty 
with this approach is that the degrees of the cofactors Vi at 
each reduction step are poorly controlled: we lack the Eu- 
clidean division step and we reduce poles at finite distance 
at the cost of making worse the pole at infinity. This diffi- 
culty is overcome by working in the projective space. The 
translation between affine and projective is discussed more 
precisely in Section 7. 

Now, assume that a and / are homogeneous polynomi- 
als in L[x] = L[xo, . . . , x n ], with / of degree d. A cen- 
tral role is played by Jac /, the Jacobian ideal of f, the 
ideal generated by the partial derivatives dof, . . . , d n f. Note 
that since / is homogeneous, Euler's relation, which asserts 
that / equals ^ yj™_ Xidif implies that / £ Jac /. 

We now decompose a as r + . Vjdjf. In contrast with 
the affine case, each nonzero v% can be chosen homogeneous 
of degree precisely dega — degdi f. If t > 1, we obtain 



F = 



1 n rj n 

*t + fi-i + / ; ' 

\ / i— 

Fi 



(£-l)f* 



(2) 



The rational function r/f is the reduced form of F and is 



If r is not zero, the order of the pole need not decrease, 
contrary to the affine case, but r is reduced to a normal 
form modulo Jac /; this will help us obtain the confinement 
property, see Proposition 2. The reduction process proceeds 
recursively on Fx, which has pole order I — 1, and stops 
when 1=1. This procedure is summarized in Algorithm 1. 



Input F = a/f a rational function in Xo, ■ ■ ■ , x n 
Output [F] such that there exist Ao, ■ ■ ■ ,A n rational 
functions such that F = [F] + diAi 

Precompute a Grobner basis G for (dof, ■ ■ ■ , d n f) 
procedure Reduce(ci//^) 
if I = 1 then return a/f 

Decompose a as r + Vidif using G 

TP, 4 L_ V 9 i Vi 

return + Reduce(Ji) 



Algorithm 1. Griffiths-Dwork reduction 

3. PROPERTIES OF THE GRIFFITHS- 
DWORK REDUCTION 

Let / in L[x] be a homogeneous polynomial of degree d, 
where L is a field of characteristic zero. It is clear that the 
reduction procedure satisfies the soundness and the linear- 
ity properties. Analogues of confinement and normalization 
hold under the following regularity hypothesis: 

L[x]/Jac/ is finite-dimensional over L. (H) 

Geometrically, this hypothesis means that the hypersurfacc 
defined by / in P™ is smooth. In particular / is irreducible. 

The ring of rational functions in L(x) whose denominator 
is a power of / is denoted L[x, 4]. Let L[x, j] p denote the 
subspace of homogeneous functions of degree p, i.e. of F 
in L[x, 4] such that -F(Ax) equals \ p F(x). Note that each 
derivation di induces a map from L[x, j] p to L[x, y] P -i. Let 
Df denote the subspace of L[x, 4] of rational functions ^\ diAi 

for some Ai in L[x, 4]-n- A major character of this study 
is the quotient space L[x, j]- n -i/Df, denoted Hf T . 

The reduced form of F in L[x, 4]_„_i is denoted [F]. It 
is by definition the output of the algorithm Reduce. It 
depends on a choice of a Grobner basis of Jac /, but its 
vanishing does not, see Theorem 1. 

The choice of the space L[x, 4]-n-i and the degree — n— 1 
may seem arbitrary. It is motivated by it being isomorphic 
to the space of regular differential n- forms on P n \ V(f). The 
evaluation of xq to 1 is the restriction map to A" \ V(f). 
The space Hf 1 is the nth de Rham cohomology space of the 
algebraic variety P n \ V(f) over k(t). 

Theorem 1 (Griffiths [12, §4]). /// satisfies Hy- 
pothesis (H), then for all F in L[x, 4]_ n -i, the reduced 
form [F] vanishes if and only if F is in Df. 

Theorem 1 gives access to the dimension of Hy . Let A be 
the finite dimensional vector space L[x]/Jac /. For a positive 
integer £, let Ai denote the linear subspace of A generated 
by homogeneous polynomials of degree id — (n + 1). Let B 
denote (BtAe- Finally, for £ > let (ge,i)i^i^n e be a basis 
of Ai, with ni = dim_L At- 

Proposition 2. Under Hypothesis (H) , the family of ra- 
tional functions (gi^/ f) Q<[ ^ induces a basis of Hf T . 

Proof. Suppose there exists a linear dependency relation 
between the guff modulo Df, that is y~] , . ui^gi^l f , de- 
noted F, lies in Df for some elements ue ti of L, not all zero. 
Let £o be the maximum £ such that uu is not zero for at least 



one i. By Theorem 1, [F] = so that y\ . Ui^gi ti f l ° 1 , the 
numerator of F, lies in Jac/. Since / itself is in Jac/, the 
sum ^\ Ui Q ^gi 0t i is in Jac /, which contradicts the fact that 
the gi 0} i are a basis of At . Thus the gi.i/f t form a free 
family. 

To prove that this family generates Hf", we first notice 
that the family of all the fractions [F], for F in L[x, 4]- n -l, 
generates Hf T since [F] equals F modulo Df. Now we 
assume for a moment that each ge,i is reduced with re- 
spect to a Grobner basis G of Jac /. Then each polyno- 
mial of L[x] of degree Id — n — 1 which is reduced with 
respect to G is a linear combination of the gz^. Thus for 
all F — a/f 1 in L[x, 4]_„_i, the reduction [F] is in the 

span of all the gi,i/f e . This makes the gi,i/f e a system of 
generators of Hf and by the previous paragraph a basis of 
it. Thus Hf T has the same dimension as B and any free fam- 
ily of Hf r of cardinal dimz, B is a basis of Hf T . In particular, 
the ge.i/ f l form a basis even if the gej are not reduced with 
respect to G. □ 

Corollary 3. Under Hypothesis (H), Hf 1 has dimen- 
sion 

PROOF. It has the dimension of B, see [19, thm. 8.4] for 
its computation. The inequality is clear. □ 

4. CREATIVE TELESCOPING 

We now introduce an algorithm, based on the Grifnths- 
Dwork reduction, that computes a telescoper of a rational 
function under Hypothesis (H). 

In Eq. (1), the telescoper T is said to have a regular certifi- 
cate if the irreducible factors of the denominators of the Ai's, 
as rational function over k(t), divide the denominator of F; 
in other words, the Ai's have no pole outside those of F, 
over k(t). Algorithm 2, described in §4.2, returns the tele- 
scoper of minimal order having regular certificate. For the 
application of creative telescoping to integration, this class 
of telescoper is more interesting than the general one; that 
is the object of §4.1. 

4.1 Telescopers with regular certificate 

Back to the afiine case, let F(t, xi, . ■ ■ , x n ) be a rational func- 
tion over C and 7 be a n-cycle in C" over which F has no 
pole for a generic t in C. A common use of creative telescop- 
ing is the computation of a differential equation satisfied by 
the one-parameter integral I(t) = J _Fdx. As mentioned in 
the introduction, it is not always possible to deduce from the 
telescoping equation (1) that T(I) vanishes. It may happen 
that the polar locus of the certificate meets 7 for all ( 6 C, 
and so J diAidn may not be zero. An example of this phe- 
nomenon is given by Picard [23] for a bivariate algebraic 
function and translated here into a rational example, using 
the method in [4, Lemma 4]: 

z' 2 - Pt{x)P t {y) = X (x-y){z 2 -k(x)Pt(y)) + 

a 2Pt{y) , Q 3(x 2 +y 2 )z 

V ( X -y)(z2-P t (x)P t (y)) Z (x-y)(z2-P t (x)P t (y)) ' K > 

where Pt (u) = u 3 +t. Note the factor x — y in the denom- 
inator of the certificate. Let F denote the left-hand side. 



Input F = a/f e a rational function in L[x, j]_ n _i, 

with / satisfying (H) 
Output T(t,dt) an operator such that T(F) = JJ\ ftAi for 

some rational functions Aj 

procedure Telesc(F) 
Go «- Reduce(F) 
i <- 
loop 

if rankt(Go, . . . , Gi) < i + 1 then 

solve 2fe=o a fcCfc = Gi w.r.t. a , . . . , tti-i in L 

return d\ - J2 k a kdt 
else 

Gi+i <— REDUCE(9 t Gi) 
i i + 1 



Algorithm 2. Creative telescoping, regular case 

In this case, the operator 1 is a telescoper of F, however 
there exists a 3-cycle 7 on which F has no pole and such 
that Fdx is not zero. It is impossible to find a regular 
certificate for the telescoper 1. 

Nevertheless, a differential equation for I(t) can be ob- 
tained in two ways. First, one can carefully study the inte- 
gral JJ\ J 9iAidx and compute a differential equation for 
it. Usually this includes the analysis of the poles of the Ai's, 
and the search of a telescoper for some rational function 
with one variable less. The second way is to find a tele- 
scoper for F such that the certificate does not contains new 
poles, a telescoper with regular certificate. Contrary to the 
telescoper (3), the operator dt is a telescoper with regular 
certificate: 

d t F = d4§F)+d y (l t F)+d z (iF). 

This proves that 8 t I = 0. More generally we have: 

Proposition 4. IfTe C(t){d t ) is a telescoper of F with 
regular certificate, then T(I) is zero. 

In this case, the certificate itself is not needed to prove the 
conclusion, its existence and regularity are sufficient. The 
Griffiths-Dwork method always produces a telescoper with 
regular certificate, see Eq. (2). 

4.2 Algorithm 

In this section L is k(t) for some field k and / is a homo- 
geneous polynomial over L of degree d satisfying Hypothe- 
sis (H). For F a rational function in L[x, j]_ n _i we want 
to find a nonzero operator T in L(dt) such that T(F) lies 
in Df. Algorithm 2 describes the procedure Telesc that 
outputs such a telescoper. 

Proposition 5. Algorithm 2 terminates and outputs the 
minimal telescoper with regular certificate of F. 

Proof. The sequence (Gfe) is defined by Go = [F] and 
the recurrence relation Gk+i = [dtGk]- We show by induc- 
tion that for all k the fraction Gk equals [dt F]. It is clear 
for k = 0. Assume that Gfe equals [8% F\. By the soundness 
of the reduction the operator Gfe — dtF is in Df. And then 
so is dtGk — dt +1 F since dt commutes with the di's. By The- 
orem 1 and linearity, this implies that [dtGk] equals [8^ 1 F]. 



At the ith step of the loop the algorithm is looking for a 
linear relation between [F], . . . , [81 F]. By Theorem I, there 
is one if and only if there is a telescoper with regular cer- 
tificate of order i. If there such a relation, it computes it 
and returns the corresponding telescoper. By the confine- 
ment property of Proposition 2, the algorithm terminates. 
By construction the telescoper admits a regular certificate, 
see Identity (2). □ 

5. EFFECTIVE BOUNDS FOR CREATIVE 
TELESCOPING 

We now review the steps of the algorithm with the aim of 
bounding the degrees and orders of all polynomials and op- 
erators that are constructed. This is then used in the next 
section to assess the complexity of this approach. 

For the needs of Section 7, we track the degrees not only 
with respect to the parameter t but also to another free vari- 
able e of the base field. In other words, we assume that L 
is k(t,s). For p a polynomial in k[t, el, the bidegree (deg t p,deg ( 
of p is denoted 8{p). If p= X^P/x 7 is a polynomial in t, e 
and x, then 5(p) denotes the supremum of the S(pi)'s, com- 
ponent by component. 

Theorem 6. Let f S L[x] be homogeneous of degree d 
satisfying (H). Let a/f e in L[x, j]_ n _i be a rational func- 
tion, with a polynomial in t and e. The minimal telescoper 
with regular certificate of a/f e has order at most d n and 
degree 

O (d n 5(a) + (£d 2n + d Zn ) e n S(f)) , 
uniformly in all the parameters. 

The first part of the theorem is a direct consequence of the 
confinement property of Corollary 3. We now study more 
precisely the decomposition used in Algorithm 1 in order to 
control the degree of the certificate and complete the proof. 

The notation a(n) = C(b(n)), for a tuple n, means that 
there exists a G > such that for all n ^ 1, with at most a 
finite number of exceptions, we have a(n) ^ Cb(n). The no- 
tation a(n) = (5(6(n)) means that o(n) = 0(b{n) log fc 6(n)) 
for some integer k. We emphasize that when there are sev- 
eral parameters in a O, the constant is uniform in all the 
parameters and there is at most a finite number of excep- 
tions. 

5.1 Reduction modulo the Jacobian ideal 

An important ingredient of the Griffiths-Dwork reduction 
is the computation of a decomposition r + y^ . Uidif of a 
homogeneous polynomial a. This can be done by means of 
a Grobner basis of Jac/, but instead of following the steps 
of a Grobner basis algorithm, we cast the computation into 
a linear algebra framework using Macaulay's matrices, for 
which Cramer's rule and Hadamard's bound can then be 
used. While not strictly equivalent, both methods ensure 
that r depends linearly on a and vanishes when a is in Jac /. 
For a positive integer q, let ip q denote the linear map 

ip q : (ui) G L[x]"±2_„ — > ^mdif G L[x] 9 _„_i. 

Let Mat y> q be the matrix of tp q in a monomial basis. It 
has dimension R q x C q , where R q denotes ( q ~ 1 ) and C q de- 
notes (n+1) ( q ~ d ) , and we note for future use that C q ^ R q 



for all positive integers n and d > 2. Up to a change of order- 
ing of the bases of the domain and codomain, Mat tp q has the 
form ( c B )> wnere ^4 is a square submatrix of maximal rank. 
Note that D is necessarily CA~ 1 B. Then, the endomor- 
phism ip q defined by the matrix (^p 1 p) satisfies ip q ip q ip q = 
ip q ; it is called a split of ip q . It depends on the choice of 
the maximal rank minor. The map id— <p q ip q , denoted 7r q , 
performs the reduction in degree q — n — 1: it is idempotent; 
if a of degree q — n — 1 is in Jac / then it equals Lp q (b) for 
some b and thus ir q (a) vanishes; and for all a in L[x] q _ n _i 
it gives a decomposition 



a = 7T«j(a) + ^>g(a) 



A/- 



Under Hypothesis (H), the map is surjective when q is 
at least (n + l)d — n. Let D denote this bound, known as 
Macaulay's bound [18, chap. 1; 16, corollaire, p. 169]. 

For q larger than D, & split of ip q can be obtained from 
a split ^jd of ifD in the following way. Let S be the set of 
monomials in x of total degree q — D. Choose a linear map /i 
from L[x] 9 _ n _i to L[x.} D _ n _ 1 such that each a in L[x] 9 _ n _i 
equals ^2 meS rn(j, m (a). Then a split of ifd is defined by 



mgS 

Let q be a positive integer and let E q be the least common 
multiple of the denominators of the entries of Mat ip q . The 
entries of Mat ip q and Mat ir q are rational functions of the 
form p/E q , with p polynomial. Let Se denote the supremum 
of all S(p) and all 8(E q ), for q € N \ {0}. 

Proposition 7. The supremum 5e is finite and bounded 
above by e n d n 8(f). Moreover, if q > D then E q equals Ed- 

Proof. Assume first that q > D. In this case, the en- 
tries of Mat?/><3 are entries of Mat?/)_D and n q is zero. Thus 
the inequalities will follow from the case where q ^ D. 
Let Mat tpq and Mat n q be written respectively as N / E q 
and P/E q with N and P polynomial matrices. Let r be the 
rank of ip q . The maximal rank minor A in the construction 
of ip q has dimension r. Cramer's rule and Hadamard's bound 
ensure that 5(N) is at most (r — l)5(f) and that 5(E q ) is at 
most r5(f). Since P equals E q id —(Mat (p q )N and <5(Mat tpd) 
equals 5(f), the degree 5(P) is also at most rS(f). 

Next, r is bounded by R q , the row dimension of Mat 0^. 
Since q ^ D, we have R q ^ Rn and we conclude using the 
following inequality, with p^nan integer: 



pe 

r+T 



□ 



Algorithm 3 is a slightly modified version of Algorithm 1 
which uses the construction above. Its output is in general 
not equal to the output of the former version, for any mono- 
mial order, but of course it satisfies Theorem 1. In particular 
the output of the algorithm Telesc does not depend on the 
reduction method in Reduce. From now on the brackets [•] 
denote the output of Algorithm 3. 

5.2 Degree bounds for the reduction 

Proposition 8. Let a/f 1 e L[x, y]_ n _i, with a a poly- 
nomial in t and e. Then 



Input F = a/f e a rational function in xo, ■ ■ ■ , x n , with / 
of degree d 

Output [F] such that there exist Ao, . . . ,A n rational 
functions such that F = [F] + diAi 

For all 1 ^ i ^ £, precompute a split ipid of tfid (§5.1) 
procedure R,EDUCE(a// £ ) 
if £ — 1 then return a/f 

1 ditpid(a)i 



i-i^ ft 



return 



+ Reduce(Fi) 



Algorithm 3. Griffiths-Dwork reduction, linear algebra variant 



where Pe = I"!— i an ^ ° k in -^[ x ]*<2— »— i is a polynomial 
in t and e such that 5(bu) 8(a) + £5e, for 1 ^ k ^ n. 

Proof. Using Algorithm 3, we obtain 



Eidf e 



1 

Eu 



f 



where g and p are polynomials in x, t and e, with S(p) 
and 5(g) at most 5(a) + 5e- Induction over I yields 



= y — 



Pi: 



, ^ i ll k E jd 



with pk polynomials such that 5(pt) ^ 5(a) + (I — k + 1)5e- 
For k > n, and hence kd > D, the map iv^d is and thus so 
is pk- Thus 



n 



E, 



y± Pk II 

k=l 



k-1 
3=1 



E 



f k 



□ 



This proposition applied to d\(aj f ) asserts that 

1 b' 



ft 



—y 



bj,k 

IF 



p ^ ^ f k p ^ f r ' 



for some polynomials bi t k and b[ such that 

6(h, k ) ^5(a)+i5(f) + (i + £)5 E , 
and 5(bi) < 5(a) + (i + n)5(f) + (i + £)c 

5.3 Degree bounds for the telescoper 



(4) 



(5) 
(6) 



Proposition 9. Let T = X^=o Ci< ^' ™^ coefficients ci 
in k[t,e], be the minimal telescoper with regular certificate 
of a/f. Then 

5(a) < r5(a) + (r 2 + rl) e n d n 5(f). 

Proof. The operator T is the output of TELESC(a// <? ). 
The rational functions Ci jc r form the unique solution to the 
following system of inhomogeneous linear equations over L, 
with the Yi's as unknown variables: 



We write each 6 ijfc in (4) as J^mGS bi,k,mm, where S is the 
set of all monomials in the variables x of degree at most nd— 
n — 1. The previous system rewrites as 

>-- 1 , , 
VmeS,Vfce{l,...,n}, y Yt b -^i = ~^i 

* — ' ri+i il+r 
i — o 

There is a set I of r indices {(fco, mo), . . . } such that the 
square system formed by the corresponding equations ad- 
mits a unique solution. We apply Cramer's rule to this sys- 
tem. Let B be the square matrix (i>i,k<,m<)i,i> f° r * an d j 
between and r — 1. Let Bi be the matrix obtained by re- 
placing the row number i of B by the vector (6 ri fc, >mj )j ■ We 
get, after simplification of the factors Pf+* by multilinearity 
of the determinant, 

r . ^P-detBi 

- = • (7) 

Consequently, for all i, the polynomial Cj divides P ^ +z det.Bi 
and thus 

r 

<5(ci) < iS E + 

With the previous bound (5) on 8(bi) we get 

5(a) < rS(a) + r ( r + 1 ) + 5 E ) + r £6 B , 
which gives the result with Proposition 7. □ 

6. COMPLEXITY 

We assume that L is the field k(t) and we evaluate the al- 
gebraic complexity of the steps of Reduce and Telesc in 
terms of number of arithmetic operations in k. All the algo- 
rithms are deterministic. For univariate polynomial compu- 
tations, we use the quasi-optimal algorithms in [10]. For sim- 
plicity, we assume that d > 2 so that several simplifications 
occur in the inequalities since C q ^ R q and d > e ~ 2.72. 

6.1 Primitives of linear algebra 

The complexity of Algorithm 2 lies in operations on matrices 
with polynomial coefficients. Let A 6 k[t] nxm have rank r 
and coefficients of degree at most d. One can compute r, 
a basis of kei A and a maximal rank minor in 0(nmr"~ 2 d) 
operarions in k [30]. A maximal rank minor can be inverted 
in complexity <D(r 3 d) [13]. In particular, a matrix B such 
that ABA = A can be computed in (D(nmr"~ 2 d + r 3 d) op- 
erations in k, or 0(n 2 rnd), using r ^ n, m and w ^ 3. 

We implicitly represent a matrix A whose entries are uni- 
variate rational functions as the least common denomina- 
tor g of the entries, together with the polynomial matrix gA. 

6.2 Precomputation 

Algorithm 2 needs the splits tpid for i from 1 to the larger 
of n + 1 and t. Following §5.1, it is enough to compute tpid 
for i between 1 and n+1, each for a cost of 0(RidC 2 d 5(f)) op- 
erations in k, and then ipid can be obtained with no further 
arithmetic operation for i > n+1. Thus the precomputation 
needs O (e 3n d 3n 5(f)) operations in k. 

6.3 Reduction 



Let p(£,5(a)) be the complexity of the linear algebra based 
variant of the algorithm Reduce with input a rational func- 
tion a/f . The procedure first computes iptdfo,). Since ipu 
is precomputed, it is only the product of a matrix of dimen- 
sions Cid by Red with the vector of coefficients of a in a 
monomial basis. The elements of the matrix have degree 
at most 5e and the elements of the vector have degree at 
most S(a). Thus the product has complexity (D(RedCed(S(a) + 
Se))- Secondly, the procedure computes r as ned(a) know- 
ing ipi d {a); this has the same complexity. Thirdly, it com- 
putes Fi, computation whose complexity is dominated by 
the first step. And lastly it computes Reduce(Fi), which 
has complexity bounded by p{l — 1, 6(a) + 8e)- Unrolling 
the recurrence leads to 

p(£, 6(a)) = d(i (Jf^ ^ (6(a) + e5 E ?j . 

6.4 Main loop 

The computation of Go has complexity p(£, 6(a)). Next, d 
has shape given by (4), and is differentiated before being 
reduced, so that the cost of the computation of d+i is at 
most p(n + 1, 5(a) + (i + 2n)8(f) + (i + £)Se). Summing up, 
the computation of Go , . . . , G r has a complexity 

p(£, 6(a)) + 6 ((ed) 2n r (5(a) + r5(f) + (r + l)5 E )) . (8) 

During the ith step, the procedure computes the rank 
of i + 1 vectors with 0(e n d n ) coefficients of degree 5(6^) and 
computes a linear dependence relation if there is one. This 
is done in complexity O (i u '~ 1 e n d n 5(b' i )') . This step is quite 
expensive and doing it for all i up to r would ruin the com- 
plexity. It is sufficient to perform this computation only 
when i is a power of 2 so that the maximal i which is used 
is smaller than 2r. When the rank of the family is not full, 
we deduce from it the exact order r and perform the compu- 
tation in that order. Indeed, the rank over L of Go, . . . , G; 
is the least of r and i. This way, finding the rank and solv- 
ing has cost 6(r u '- 1 e n d n (6(b' r ) + 6(b' 2r ))). In view of (6) 
and since r ^ d n and w ^ 3, the complexity of that step 
is bounded by (8). Adding the cost of the precomputation 
and using the bounds of the previous section leads to the 
following complexity estimate. 

Theorem 10. Under Hypothesis (H) and assuming that 
d > 2, Algorithm Telesc run with input a/f e takes 

O ((d 5n + d 4n £ + d' in f (I) 2 ™) e 3n 6j 

arithmetic operations in k, where 5 is the larger of 5(a) 
and 5(f), uniformly in all the parameters. Asymptotically 
with I and n fixed, this is O (d 5n 5) . 

Note that while this may seem a huge complexity, it is not so 
bad when compared to the size of the output, which seems to 
be, empirically, comparable to d 3n 5, with n fixed and I = 1. 
Note also that for n = 1, the complexity improves over that 
of the algorithm based on Hermite's reduction studied in [3], 
thanks to our avoiding too many rank computations. 

7. AFFINE SINGULAR CASE 

Let L denote the field k(t). Let F a g be a rational func- 
tion in L(xi, . . . ,x„), written as a/f^g. We do not assume 
that Fan is homogeneous, nor that / a ff satisfies a regularity 
property. Let d a ff be the total degree of / a a w.r.t. x. 
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order of telesc. 
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0.4s 
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0.6s 
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140s 
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100 (204) 


0.9s 


519 (2673) 


270s 


1704 (16794) 


13h 
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Table 1. Empirical order and degree of the minimal telescoper with regular certificate of a random rational function a// 2 
in Q(t, xq, x\ , X2), with / and a homogeneous in x satisfying deg x a + 3 = 2 deg x / and 8(a) and 8(f) equal to 8; 
together with a proved upper bound (with a version without simplification of Theorem 9) and mean computation time 
(CPU time). 



In this section we show a deformation technique that reg- 
ularizes singular cases. In particular, it makes it possible to 
transfer the previous results to the general case and obtain 
the following bounds. 

Theorem 11. The function Fag admits a telescoper, with 
regular certificate, of order at most d" and degree 

0(dl I 5(a)+dl n I e n S(U)), 

where d pr is max(d a g, deg x a + n + 1). This telescoper can 
be computed in complexity 0{e Zn d^)5^, with 8 the larger 
of 5(a) and <5(/ a s)- 

The algorithm is again based on linear algebra. It is not 
difficult to see that the bit complexity too is polynomial 
in dp r . 

7.1 Homogenization and deformation 

The regularization proceeds in two steps. First, let F pr be 
the homogenization of F^fs in degree — n — 1, that is 

which we write 6// pr for some homogeneous polynomials b 
and / pr . Let d pl denote the degree of / pr ; it is given by 
Theorem 11. The degrees of b and / pr satisfy the hypothesis 
of Theorem 6, by construction, but in general / pr does not 
satisfy Hypothesis (H). (Although it does generically, as 
long as d pr equals d a ff-) We consider a new indeterminate e, 
the polynomial / rcg defined by 



n 




i=0 



and the rational function F re g defined by b/f Teg . 

Lemma 12. The polynomial f lag satisfies Hypothesis (H) 
over L(e), that is L(e)[x]/Jac / rog has finite dimension. 

Proof. It is true for e = 00, so it is generically true. □ 

Now, Theorem 6 gives bounds on the order and degree of a 
telescoper of F Teg , which is in L(e)[x, - — ]_ n _i. The end of 
the proof of Theorem 11 follows from the following. 

Proposition 13. IfT in L[e](d t ) is a telescoper of F Icg 
with regular certificate, then so is Tj £= o for F&g. 

Proof. By assumption, T(F ICS ) equals y™_ p djgj / ff^, 
for some integer p and polynomials gi in L(e) [x] . Each gi / ff e& 
can be expanded in Laurent series in e as 5Zj>jv hijS 1 f° r 
some possibly negative integer N and rational functions hij 
in L[x, yL]_„. Similarly, we can write the operator T(F Teg ) 



as Tj e= o(-Ppr) + £^2j >0 bj£ :> for some rational functions bj 

in Z/[x, — ]. Since the derivations di commute with e, it is 

clear that Tj e=0 (.Fp r ) equals ^™=o ^^o- Next, in this equal- 
ity, xo can be evaluated to 1 to give 

n 

T\ e= o(F&ft) = (dohoo)\x =i + ^ dj(hio\ xo =i). 

1=1 

Euler's relation for hoo gives (with the index 00 dropped) 

71 

(doh)\ xo =i = — dj(xih\x =i), 

i=l 

proving that (doh)\ xg =i is in Df aB . Thus, so is T| e=0 (-F a ff) 
and the proof is complete. □ 

Nevertheless, a telescoper obtained this way does not need 
to be minimal, even starting from a minimal one for the per- 
turbed function -F rog . This is unfortunate because in pres- 
ence of singularities the dimension of H^ T can collapse when 
compared to the generic order given by Corollary 3. Note 
also that in the proof of the previous Proposition, the inte- 
ger N can be bounded in d 0( - n ^ a which implies a bound on 
the size of the certificate in d 0< - n ' similar to Lipshitz's [17]. 

7.2 Algorithm and complexity 

The algorithm is based on Proposition 13. We use an evaluation- 
interpolation scheme to control the complexity. Let the op- 
erator T in k(t, e){dt) be the minimal telescoper of F Iog , writ- 
ten as dl + J2h=o T~®t- It is tne output of Telesc applied 
to F leg . We aim at computing (£ a T)| e = ! where a is such 
that this evaluation is finite and not zero. 

Proposition 13, slightly adapted, shows that T\ e=u is a 
telescoper with regular certificate of F^ g u whenever c r (t, it) 
is not zero, even if fle g u does not satisfy (H). When it does, 
the specialization gives the minimal one: 

Lemma 14. If flt g u satisfies hypothesis (H) and if u does 
not cancel c r , then Ti s=tI is the minimal telescoper with reg- 
ular certificate of F}^ g u . 

Proof. We use the notation of Section 5, replacing / 
by / rcg and L by L(e). The operator T is the output of 
Algorithm 2 applied to F rcg . Since fll g u satisfies (H), for 
all d the matrix Mati/>d, with coefficients in L[e], has the 
same rank as its specialization with e = u [18, §58]. Thus, 
to compute the splits ibd we can choose maximal rank mi- 
nors of Mat ifid that are also maximal rank minors of the 
specialization. When doing so, the reduction [•] commutes 
with the evaluation ■| £=u . In particular, the polynomials E q 
do not vanish for e = u. 



In the proof of Prop. 9, Eq. (7) shows that c r , the leading 
coefficient of T, divides Pi +r AetB. The polynomial Pe+ r is 
a product of several Ed's, in particular Pt.+ r \e=u is not zero. 
Since c r | e=u 7^ 0, the determinant of B\ e=u is not zero either. 
Looking at the definition of B in the proof of Proposition 9, 
this implies that the [dlF les ]i s=u , for i between and r — 1 
are free over L(s). In particular, a telescoper with regular 
certificate of F}^ u has order at least r. Since T| e = u is a 
telescoper of order is r, it is the minimal one. □ 

We now present the algorithm. Let N be e" (d 3n + d 2n + 
d n ). By Proposition 9, the polynomials Ck have degree 
at most in e, and at most N5 in t. Choose U a set 
of 4JV + 1 elements of k. Determine U' the set of elements u 
of U such that f\lg U satisfies (H). This step has complex- 
ity 6((ed) n "8\U\): The polynomial f\%= u satisfies (H) if 
and only if (Mat ipo)\e=u is full rank. In particular, if f\%g W 
does not satisfy (H), then Er>\ e=u vanishes. The polyno- 
mial Ed has degree at most e"d n in e, by Proposition 7, 
so U \ U' has at most e n d n elements. For each u in U' , 
compute TELESC(/Jcg = ") with leading coefficient normalized 
to 1, denoted T u . This step has complexity <D(d 5n e 3n S\U'\), 
by Theorem 10. Determine the subset U" of U' where the or- 
der of T u is maximal. By Lemma 14, the complement U'\U" 
is formed by u such that c r (t,u) = 0. It has at most N ele- 
ments since c r has degree at most N ins. For all u in U" the 
operators T u and Tj e=u coincide. Thus U" has at most 2A r +l 
elements. 

The r rational functions ° fc (t'o) can ' 3e com P u ted using 
Lemma 15 in total complexity <D(N 2 rS). If c r (i,0) is zero, 
we look for the positive integer a such that the functions e a ^4 
are finite for e = but not zero for at least one k. The in- 
teger a is at most N and thus can be found with a binary 
search, using at most log 2 N + 1 times Lemma 15. 

Lemma 15. Let R in k(x,y) be written P/Q, with P and Q 
polynomials of degree less than d x in x and d y in y. Given 
evaluations R(x,v), for 2d y + 1 elements v of k, the func- 
tion R(x,0) (or 00 if Q(x,Q) vanishes) can be computed us- 
ing 0(d x d y ) arithmetic operations in k. 

Proof. Let V the set of evaluation points. Choose a 
set U of 2d x + 1 points of k. Compute R(u, v) for u 6 U 
and v £ V in (D(d x d y ) operations. Note that there is no need 
to check that the elements of U are not poles of the R(x, v): 
univariate rational reconstruction can handle that. Use uni- 
variate rational reconstruction to compute R(u, y) , for u 6 
U, in complexity 0(d y \U\) operations. Reconstruct R(x,0) 
in complexity <D(d x ) from the evaluations R(u, 0). □ 

8. EXPERIMENTS 

A basic implementation of the algorithm Telesc has been 
written in Maple 16. As it uses only Maple primitives to 
compute with polynomial matrices, it is certainly too basic 
to reflect the complexity given in Theorem 10. 

Table 1 presents empirical results for some generic rational 
functions, with n = 2. The bound on the order are generi- 
cally exact as expected; however the bound on the degree are 
not very sharp. For n — 1 and 5(a) fixed, a careful study [3] 
proves that the degree of the minimal telescoper is 0(d 2 5), 
which is smaller than the 0(d i 8) given by Theorem 6. Anal- 
ogy, as well as numerical evidence and theoretical clues, lead 



us to think that for general n, the asymptotic behavior can 
be improved from 0(d in S) to 0(d 2n S). 

The relative cost of each step of Algorithm 2 in the com- 
putation of telescopers of Table 1, on the example of the 
telescoper of degree 12 and degree 1092 of a generic func- 
tion a/f 2 as described in Table 1, that is computed in about 7 
hours breaks down as follows: The computation of splits of 
Macaulay matrices takes about 1% of the time, the reduction 
steps about 40%, and the final solving about 60% of the time. 
More efficient matrix multiplication and system resolution 
over univariate polynomials could improve speed dramati- 
cally. We have not been able to compute more than the first 
column of Table 1 with methods and programs in [15,4]. 

On the other hand, the regularity hypothesis (H) is re- 
strictive in applications: Even though generic polynomials 
satisfy this hypothesis, examples with physical or combinato- 
rial meaning usually do not. The method shown in Section 7 
is only of a theoretical interest. By contrast, the algorithm 
for the regular case is very efficient in practice. 
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